/*
 * Main.cpp
 *
 *  Created on: 13 Jul 2011
 *      Author: allan
 */

#include <iostream>

#include "System.h"
//#include "SystemBuilder.h"
#include "EquilibriumSolver.h"
#include "MultirootSolver.h"

int main()
{
//	double T = 25.0, P = 100.0;
//	
//	StringVector knownSpecies  = {"H2O", "H[+]", "Na[+]", "Cl[-]", "NaCl(aq)", "OH[-]"};
//	
//	DoubleVector n = {55.508, 1.0E-7, 1.0, 1.0, 0.01, 1.0E-7};
//	
//	const System system = SystemBuilder::GenerateSystemNoSpeciesGeneration(knownSpecies);
//	
//	EquilibriumSolver speciationSolver(system, {"NaCl(aq)", "OH[-]"});
//	
//	speciationSolver.SpecifyChargeBalance();
//	speciationSolver.SpecifySpeciesMoles("H2O", 55.508);
//	speciationSolver.SpecifySpeciesMoles("Na[+]", 1.0);
//	speciationSolver.SpecifySpeciesMoles("Cl[-]", 1.0);
//	speciationSolver.SetSolverOptions(1.0E-6, 50);
//	
//	n = speciationSolver.Solve(T, P, n);
//	
//	cout << setprecision(15);
//	
//	cout << n << endl;
//	
//	DoubleVector a = system.Activities(T, P, n);
//	
//	cout << log10(system.GetReactions()[4].EquilibriumConstant(T, P)) << endl;
//	cout << log10(system.GetReactions()[5].EquilibriumConstant(T, P)) << endl;
//	
//	cout << log10(system.GetReactions()[4].ReactionQuotient(a)) << endl;
//	cout << log10(system.GetReactions()[5].ReactionQuotient(a)) << endl;
//	
//	cout << log10(system.GetReactions()[4].EquilibriumConstant(T, P)/system.GetReactions()[4].ReactionQuotient(a)) << endl;
//	cout << log10(system.GetReactions()[5].EquilibriumConstant(T, P)/system.GetReactions()[5].ReactionQuotient(a)) << endl;
//	
	return 0;
}

//void func(const VectorXd& x, VectorXd& F)
//{
//	F[0] = x[0]*x[0] + x[1]*x[1] - 34.0;
//	F[1] = x[0]*x[1] - 15.0;
//}
//
//void jac(const VectorXd& x, MatrixXd& J)
//{
//	J(0,0) = 2*x[0];
//	J(0,1) = 2*x[1];
//	J(1,0) = x[1];
//	J(1,1) = x[0];
//}
//
//void funcjac(const VectorXd& x, VectorXd& F, MatrixXd& J)
//{
//	func(x, F);
//	jac(x, J);
//}
//
//int main()
//{
//	MultirootSolver solver;
//	
//	solver.SetSolver(2, func, jac, funcjac);
//	
//	VectorXd x(2); x << 2, 8;
//	
//	solver.Solve(x);
//	
//	cout << x;
//	
////	cout << "Allan" << endl;
//	
//	return 0;
//}


//	double T = 25.0, P = 100.0;
//	
//	StringVector knownSpecies  = {"H2O", "H[+]", "Na[+]", "Cl[-]", "NaCl(aq)", "OH[-]"};
//	
//	Phase aqueousPhase = Phase::CreateAqueousPhase(knownSpecies);
//	
//	DoubleVector n = {55.508, 1.0E-7, 1.0, 1.0, 0.01, 1.0E-7};
//	
//	DoubleVector a = aqueousPhase.Activities(T, P, n);
//	
//	cout << a << endl;
//	
//	const System system = SystemBuilder::GenerateSystemNoSpeciesGeneration(knownSpecies);
//	//cout << system << endl;
//	
//	cout << system.Activities(T, P, n) << endl;
//	
//	BOOST_FOREACH(const Reaction& reaction, system.GetReactions())
//	{
//		cout << log10(reaction.EquilibriumConstant(T, P)) << endl;
//		cout << log10(reaction.ReactionQuotient(a)) << endl;
//		cout << endl;
//	}
